// panda.cpp : Defines the entry point for the console application.
//

#if _MSC_VER > 1000
#pragma warning( disable :  4511 4512 4786)
// nasty hack for VC++ non-compliance with standard C++ for loop declaration scope
#define for if (0); else for
#endif

#include <stdio.h>
#include <math.h>
#include <string>
#include <set>
#include <iostream>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <iterator>
#include "boost/tuple/tuple.hpp"
#define _MY_STL_MACROS
#include "components/stlhelper.h"
#include "components/perlutil.h"
#include "components/popen/popen.h"
#include "components/progress_indicator.h"
#include <assert.h>
#include "blosum_data.h"
#include "components/portable_timer.h"
#include "get_arg.h"

using std::string;
using std::cerr;
using std::cin;
using std::vector;
using std::fill;
using std::cout;
using std::ifstream;
using std::ofstream;
using std::ostream;

#define BLOSUM80
#ifdef BLOSUM60
#define GAP_EXIST_PENALTY -11
#define BLOSUM_MATRIX blosum62
#else
#ifdef BLOSUM80
#define GAP_EXIST_PENALTY -10
#define BLOSUM_MATRIX blosum80
#endif
#endif
// global variables!!
// This must be a perl programme in disguise!!
namespace
{
	string my_path;
	string tmp_file_path;
	string yn00_exe_path;
	string paml_path;
	char lookup_table[256][256];
	char gap_penalty[256][256];
}

string GetPathPart(const string& path)
{
	string::size_type p1, p2, p;
	p1 = path.rfind('/');
	p2 = path.rfind('\\');

	// neither found : return this directory
	if (p2 == p1)
		return ".";

	if (p1 != string::npos)
	{
		if (p2 != string::npos)
			p = std::max(p1, p2);
		else
			p = p1;
	}
	else
		p = p2;
	return path.substr(0, p);
}













//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//	loci code
//
struct t_segment
{
	unsigned beg;
	unsigned end;
	t_segment(unsigned beg_, unsigned end_)
		: beg(beg_), end(end_){}
	t_segment(const string& data)
	{
		beg = atoi(data.c_str()) - 1;
		end = atoi(data.c_str() + data.find(':') + 1) - 1;
		if (beg  > end)
			beg = end = 0;
	}
};
typedef vector<t_segment> t_loci;



//________________________________________________________________________________________

//	merge_overlapping_loci

//________________________________________________________________________________________
void merge_overlapping_loci(t_loci& loci,
							t_loci& complement_loci,
							unsigned aln_beg, unsigned aln_end)
{
	vector<char> bits(aln_end, '\x0');
	typedef vector<char>::iterator iter_pos;

	// mark bits of all masked regions
	for (unsigned i = 0; i < loci.size(); ++i)
	{
		fill(	bits.begin() + std::min(loci[i].beg, aln_end),
				bits.begin() + std::min(loci[i].end + 1, aln_end),
				0x1);
	}

	// non-overlapping masked loci
	t_loci loci2;

	// get masked regions
	iter_pos first = find(bits.begin() + aln_beg, bits.end(), 0x1);
	while (first != bits.end())
	{
		iter_pos last = find(first+1, bits.end(), 0x0);
		loci2.push_back(t_segment(first - bits.begin(),
								   	   last  - bits.begin() - 1));
		first = find(last, bits.end(), 0x1);
	}
	loci2.swap(loci);

	// get non-masked regions
	first = find(bits.begin() + aln_beg, bits.end(), 0x0);
	while (first != bits.end())
	{
		iter_pos last = find(first+1, bits.end(), 0x1);
		complement_loci.push_back(t_segment(first - bits.begin(),
												 last  - bits.begin() - 1));
		first = find(last, bits.end(), 0x0);
	}
}

//________________________________________________________________________________________

//	subtract_overlapping_loci

//________________________________________________________________________________________
void subtract_loci(	const t_loci& original_loci, const t_loci& remove_loci,
					t_loci& resulting_loci,
			   		unsigned aln_beg, unsigned aln_end)
{
	resulting_loci.clear();
	vector<char> bits(aln_end, '\x0');
	typedef vector<char>::iterator iter_pos;

	// mark bits of all loci
	for (unsigned i = 0; i < original_loci.size(); ++i)
	{
		fill(	bits.begin() + std::min(original_loci[i].beg, aln_end),
				bits.begin() + std::min(original_loci[i].end + 1, aln_end),
				0x1);
	}

	// remove bits of all loci
	for (unsigned i = 0; i < remove_loci.size(); ++i)
	{
		fill(	bits.begin() + std::min(remove_loci[i].beg, aln_end),
				bits.begin() + std::min(remove_loci[i].end + 1, aln_end),
				0x0);
	}

	// get resulting loci
	iter_pos first = find(bits.begin() + aln_beg, bits.end(), 0x1);
	while (first != bits.end())
	{
		iter_pos last = find(first+1, bits.end(), 0x0);
		resulting_loci.push_back(t_segment(first - bits.begin(),
			   						   last  - bits.begin() - 1));
		first = find(last, bits.end(), 0x1);
	}
}


//________________________________________________________________________________________

//	 extract_seq_pair_with_loci

//________________________________________________________________________________________
inline string extract_seq_with_loci(const std::string& seq, const t_loci& loci)
{
	string nrvo;
	for (unsigned i=0; i < loci.size(); ++i)
		nrvo += seq.substr(loci[i].beg,  loci[i].end + 1 - loci[i].beg);
	return nrvo;
}
inline void extract_seq_pair_with_loci(vector<string>& sequences,
										const std::string& seq1,
										const std::string& seq2,
										const t_loci& loci)
{
	// add sequence for mouse
	string seq_region;
	for (unsigned i=0; i < loci.size(); ++i)
		seq_region += seq1.substr(loci[i].beg,  loci[i].end + 1 - loci[i].beg);
	sequences.push_back(seq_region);

	// add sequence for human
	seq_region.clear();
	for (unsigned i=0; i < loci.size(); ++i)
		seq_region += seq2.substr(loci[i].beg,  loci[i].end + 1 - loci[i].beg);
	sequences.push_back(seq_region);
}

//________________________________________________________________________________________

//	 extract_cdna_pair_with_loci

//________________________________________________________________________________________
struct op_prot_loci_to_dna
{
	void operator() (t_segment& a)
	{
		a.beg *= 3;
		a.end = (a.end + 1) * 3 - 1;
	}
};

inline void extract_cdna_pair_with_loci(vector<string>& sequences,
										const std::string& seq1,
										const std::string& seq2,
										const t_loci& prot_loci)
{
	// prot loci to cdna loci
	t_loci loci = prot_loci;
	for_each (beg_to_end(loci), op_prot_loci_to_dna());

	// share processing with protein function
	extract_seq_pair_with_loci(sequences, seq1, seq2, loci);
}

//
//	loci code
//
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888






//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	bitscore

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//________________________________________________________________________________________

//	get_seq_bitscore

//		calculate blosum raw bitscores for one stretch of a sequence pair

//________________________________________________________________________________________
// helper struct
struct op_calc_bitscore
{
	char	(*gap_penalty)[256];
	char	(*lookup_table)[256];
	int	bitscore;
	char prev_a;
	char prev_b;
	op_calc_bitscore(char gap_penalty_[256][256], char lookup_table_[256][256]):
			gap_penalty(gap_penalty_), lookup_table(lookup_table_),
			bitscore(0), prev_a('-'), prev_b('-'){}
	void operator()(char a, char b)
	{
		// lookup bitscore
		bitscore += lookup_table[a][b];

		// gap opening penalty
		bitscore += gap_penalty[prev_a][a];
		bitscore += gap_penalty[prev_b][b];
		// save current char to calculate gap opening penalty
		prev_a = a;
		prev_b = b;

	}
};

int get_seq_bitscore(const std::string& seq1, const std::string& seq2,
							unsigned beg, unsigned end, unsigned& len)
{
	if (!seq1.length())
		return 0;
	// skip past initial '-'
	beg = seq1.find_first_not_of('-', beg);
	beg = seq2.find_first_not_of('-', beg);

	// skip past final '-'
	end = seq1.find_last_not_of('-',  end);
	end = seq2.find_last_not_of('-',  end);
	if (beg == string::npos ||
		end == string::npos || end < beg)
		return 0;
	assert(beg >= 0);
	assert(end < seq1.length());
	assert(seq2.length() == seq1.length());
	len += end + 1 - beg;
	return (for_each_pair(	seq1.begin() + beg, seq1.begin() + end + 1,
							seq2.begin() + beg,
							op_calc_bitscore(gap_penalty, lookup_table))
			).bitscore;
}

//________________________________________________________________________________________

//	get_all_regions_bitscores

//________________________________________________________________________________________
// helper function
void print_bitscore_over_loci(const std::string& seq1, const std::string& seq2,
								const t_loci& loci)
{
	double bitscore = 0.0;
	unsigned len = 0;
	for (unsigned i=0; i < loci.size(); ++i)
		bitscore += get_seq_bitscore(seq1, seq2, loci[i].beg, loci[i].end, len);
	if (len)
		cout << "\t" << (bitscore / len);
	else
		cout << "\t\\N";
}

void get_all_regions_bitscores(	const string& prot_m, const string& prot_h,
								const vector<t_loci>& vec_loci)

{
	// bitscore for each region
	for (unsigned i = 0; i < vec_loci.size(); ++i)
		print_bitscore_over_loci(prot_m, prot_h, vec_loci[i]);
}





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	percentage_id

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


//________________________________________________________________________________________

//	get_all_regions_percentage_id

//________________________________________________________________________________________
//	helper functions
struct t_count_identical
{
	explicit t_count_identical(unsigned i = 0): count(i){}
	unsigned count;
	void operator() (char a, char b)
	{
		if (a == b)
			++count;
	}
};

void get_all_regions_percentage_id(	vector<string>& proteins)
{
	// percentage id = -1.0 for empty strings
	for (unsigned i = 0; i < proteins.size(); i += 2)
	{
		if (proteins[i].length())
		{
			cout << "\t" << (
							 for_each_pair(	beg_to_end(proteins[i]),
									proteins[i+1].begin(),
									t_count_identical()).count
									*100.0 / proteins[i].length());
		}
		else
			cout << "\t\\N";
	}
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	KaKs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//	helper function for below
// returns whether isnan / isinf /is 0.0
inline bool print_valid_double_pg_field(double num, ostream& os)
{
	if (isnan(num) || isinf(num) || num > 50)
	{
		os << "\t\\N";
		return false;
	}
	else
	{
		os << "\t" << num;
		return true;
	}
}

//________________________________________________________________________________________

//	print_sequence_ka_ks

//________________________________________________________________________________________
bool print_sequence_ka_ks(const string& cdna_m, const string& cdna_h,
						ofstream& err_strm)
{
	if (!cdna_m.length())
	{
		cout << "\t\\N\t\\N\t\\N";
		return true;
	}
	// remove any existing results file before running analysis
	remove((my_path + "/kaks.res").c_str());

	ofstream seq_strm (tmp_file_path.c_str());
	if (!seq_strm)
	{
		cerr << "\n\nError!!!\n\tCould not open " + tmp_file_path + "\n";
		return false;
	}

	seq_strm << "2 " << cdna_m.length() << "\n";
	seq_strm << "a  \n" << cdna_m << "\n";
	seq_strm << "b  \n" << cdna_h << "\n";
	seq_strm.close();

	{
		Popen Run_yn00;
		FILE *FH_YN00_ERR;
		FILE* FH_YN00 = Run_yn00.open_stderr(yn00_exe_path.c_str(), &FH_YN00_ERR, "r");

		if (!FH_YN00)
		{
			cerr << "\n\nError!!!\n\tCould not run " << yn00_exe_path <<"\n";
			return false;
		}

		// log errors
		char err_buf[1024];
		while (size_t bytes = fread(err_buf, 1, sizeof (err_buf) , FH_YN00_ERR))
			err_strm.write(err_buf,  bytes);
		Run_yn00.close_stderr(FH_YN00, FH_YN00_ERR);
	}

	// get results
	ifstream f_results ((my_path + "/kaks.res").c_str());
	if (!f_results)
	{
		cout	<< "\t\\N\t\\N\t\\N";
		return true;
	}
	double KsKa[2];
	f_results.read((char*)KsKa, sizeof(double) * 2);
	if (f_results)
	{
		bool valid1 = print_valid_double_pg_field(KsKa[1], cout);
		bool valid2 = print_valid_double_pg_field(KsKa[0], cout);
		if (valid1 && valid2 && KsKa[0] != 0.0 && KsKa[0] != -0.0)
			print_valid_double_pg_field(KsKa[1]/KsKa[0], cout);
		else
			cout << "\t\\N";
	}
	else
		cout	<< "\t\\N\t\\N\t\\N";
	return true;
}



//________________________________________________________________________________________

//	get_all_regions_ka_ks

//________________________________________________________________________________________
bool get_all_regions_ka_ks(	const vector<string>& cDNAs)
{
	ofstream err_strm ((my_path+ "/yn00.err").c_str());
	if (!err_strm)
	{
		cerr << "\n\nError\n\tUnable to open " << my_path << "/yn00.err\n";
		return false;
	}


	for (unsigned i = 0; i < cDNAs.size(); i += 2)
	{
		if (!print_sequence_ka_ks(	cDNAs[i],
									cDNAs[i+1],
									err_strm))
			return false;

	}
	return true;
}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	matrix

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//________________________________________________________________________________________

//	get_all_regions_matrix

//________________________________________________________________________________________
struct op_calculate_matrix
{
	op_calculate_matrix(short freq_lookup_table_[256][256],	short* freq_matrix_):
			freq_lookup_table(freq_lookup_table_), freq_matrix(freq_matrix_){}
	short	(*freq_lookup_table)[256];
	short* 	freq_matrix;
	void operator() (char a, char b)
	{
		freq_matrix[freq_lookup_table[a][b]]++;
	}
};

void get_all_regions_matrix(	vector<string>& sequences,
							  	short freq_lookup_table[256][256])
{
	short freq_matrix[401];

	for (unsigned i = 0; i < sequences.size(); i += 2)
	{
		memset(freq_matrix, 0, sizeof (freq_matrix));
		for_each_pair(	beg_to_end(sequences[i]), sequences[i+1].begin(),
						op_calculate_matrix(freq_lookup_table, freq_matrix));

		//	only print out differences in the two sequences
		cout << "\t";
		for (unsigned j = 0; j < 400; ++j)
			if (freq_matrix[j] != 0 && (j / 20 != j % 20))
				cout << valid_aa_chars[j / 20] << valid_aa_chars[j % 20]
					 << freq_matrix[j] << ",";

		//	only print out identical
		cout << "\t";
		for (unsigned j = 0; j < 20; ++j)
			if (freq_matrix[j * 20 + j] != 0)
				cout << valid_aa_chars[j]
					 << freq_matrix[j * 20 + j] << ",";
	}
}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	lookup tables

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//________________________________________________________________________________________

//	initialize_freq_lookup_table

//________________________________________________________________________________________
void initialize_freq_lookup_table(short table[256][256])
{
	for (unsigned i = 0; i < 256; ++i)
		for (unsigned j = 0; j < 256; ++j)
			table [i][j] = 400;

	for (unsigned i = 0; i < 20; ++i)
		for (unsigned j = 0; j < 20; ++j)
			table [valid_aa_chars[i]][valid_aa_chars[j]] = i * 20 + j;

}

//________________________________________________________________________________________

//	initialize_bitscore_lookup_table

//________________________________________________________________________________________
void initialize_bitscore_lookup_table(char lookup_table[256][256], char gap_penalty[256][256])
{
	// bitscore lookup table
	memset (lookup_table, 0, 256 * 256 * sizeof (char));
	for (unsigned i = 0; i < sizeof (blosum_chars)-1; ++i)
		for (unsigned j = 0; j < sizeof (blosum_chars)-1; ++j)
		lookup_table[blosum_chars[i]][blosum_chars[j]] = BLOSUM_MATRIX[i][j];

	// gap extension penalty
	for (unsigned i = 0; i < sizeof(blosum_chars)-1; ++i)
	{
		lookup_table[blosum_chars[i]]['-'] = -1;
		lookup_table['-'][blosum_chars[i]] = -1;
	}

	// gap initiation penalty
	memset (gap_penalty, 0, 256 * 256 * sizeof (char));
	for (unsigned i = 0; i < sizeof(blosum_chars)-1; ++i)
		gap_penalty[blosum_chars[i]]['-'] = -11;
}

















//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	orthologs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888




struct ortholog_region
{
	std::string				gene_id1;
	std::string				gene_id2;
	std::vector<t_segment>	loci_masked;
	std::string				description;
};

void get_gene_ids(std::vector<ortholog_region>& orthologs)
{
	std::string line;
	t_progress_indicator dots(cerr, 100);

	while (cin && getline(cin, line))
	{
		std::vector<std::string> data;
		// get maximum of four parameters
        perl_split (line, "\t", data, false, 4);
		if (data.size() < 4)
		{
			cerr << "\nData not four parameters: " << line << "[" << data.size() << "]\n";
			continue;
		}
		orthologs.push_back(ortholog_region());

		// save gene ids and description
		orthologs.back().gene_id1		= data[0];
		orthologs.back().gene_id2		= data[1];
		orthologs.back().description	= data[3];
		if (! data[2].length())
			continue;
		// split on commas and put into region
		std::vector<std::string> loci_str;
		perl_split(data[2],  "|",  loci_str);
		for (unsigned i = 0; i < loci_str.size(); ++i)
			orthologs.back().loci_masked.push_back(t_segment(loci_str[i]));


		++dots;
	}
}

bool expand_gaps(const string& gaps_str, string& prot, string& cdna, const string& gene_id)
{
	// empty string encodes no gaps!
	if (! gaps_str.length())
		return true;

	// split string of gaps
	vector<string> gaps;
	perl_split(gaps_str,  "|",  gaps);
	for (int i = gaps.size()-1; i >= 0; --i)
	{
		// gaps should not come after a protein or be humungously long
		unsigned beg = atoi(gaps[i].c_str());
		unsigned len = atoi(gaps[i].c_str() + gaps[i].find(':') + 1);
		if (beg > prot.length() || len > 2000)
		{
			return false;
		}
		prot.insert(beg, string(len,  '-'));
		cdna.insert(beg * 3, string(len * 3,  '-'));
	}
	return true;
}

// eREGION {eMASKED = 1, eINVERSE = 2, eWHOLE = 4};
// enum eMASK {eMASKED, eUNMASKED, eBOTH};
bool analyze_each_ortholog(std::vector<ortholog_region>& orthologs, const string& sql_str,
							unsigned mask_region, short freq_lookup_table[256][256],
							t_pg_cursor* pg_db)
{
	ofstream err_strm ((my_path+ "/lengths.err").c_str());
	if (!err_strm)
	{
		cerr << "\n\nError\n\tUnable to open " << my_path << "/lengths.err\n";
		return false;
	}




	t_progress_indicator dots(cerr, 10);
	for (unsigned i = 0; i < orthologs.size(); ++i, ++dots)
	{
		string exec_sql = sql_str + " WHERE gene_id1 = '" + orthologs[i].gene_id1;
		exec_sql += "' AND gene_id2 = '" + orthologs[i].gene_id2 + "' LIMIT 1";

		if (!pg_db->ExecTuplesOk(exec_sql.c_str()))
		{
			err_strm <<  "\nError:\n[" << exec_sql << "]\n" << pg_db->ErrorMessage();
			return false;
		}

		if (!pg_db->Tuples())
			continue;

		string prot1 = pg_db->getvalue(0, 0);
		string prot2 = pg_db->getvalue(0, 1);
		string cdna1 = pg_db->getvalue(0, 2);
		string cdna2 = pg_db->getvalue(0, 3);
		string gaps1_str = pg_db->getvalue(0, 4);
		string gaps2_str = pg_db->getvalue(0, 5);
		unsigned align_beg = atoi(pg_db->getvalue(0, 6)) - 1;
		unsigned align_end = atoi(pg_db->getvalue(0, 7)) + align_beg;


		if (!expand_gaps(gaps1_str, prot1, cdna1, orthologs[i].gene_id1) ||
			!expand_gaps(gaps2_str, prot2, cdna2, orthologs[i].gene_id2))
		{
			//length problem dump diagnostic
			err_strm<< "\n\nError!!!\n\tGap Lengths wrong. "
						"This ortholog will be ignored.\n\t[ "
					<< orthologs[i].gene_id1	<< "/"
					<< orthologs[i].gene_id2 << "]\n gap strings: "
					<< gaps1_str << "\t" << gaps2_str << "\nProtein lengths: "
					<< prot1.length() << "\t" << prot2.length() << "\n"
					<< prot1					<< "\n"
					<< prot2					<< "\n";
			continue;
		}

		//
		// insert gaps into protein and cDNA
		//
		if (align_beg >= align_end || align_end > prot1.length())
		{
			err_strm<< "\n\nError!!!\n\tAlignment coordinates invalid. "
						"This ortholog will be ignored.\n\t[ "
					<< align_beg << "," << align_end << "]\n\n["
					<< orthologs[i].gene_id1	<< "/"
					<< orthologs[i].gene_id2 << "]\n\n";
			continue;
		}

		if (prot1.length() != prot2.length() ||
			cdna1.length() != cdna2.length() ||
			cdna1.length() != prot1.length() * 3)
		{
			err_strm<< "\n\nError!!!\n\tSequence lengths do not match. "
						"This ortholog will be ignored.\n\t[ "
					<< orthologs[i].gene_id1	<< "/"
					<< orthologs[i].gene_id2 << "]\n"
					<< prot1					<< "\n"
					<< prot2					<< "\n";
			continue;
		}

		//
		// merge overlapping regions
		//
		std::vector<t_loci> vec_loci(3);
		vec_loci[0].swap(orthologs[i].loci_masked);
		merge_overlapping_loci(vec_loci[0], vec_loci[1], align_beg, align_end);
		vec_loci[2].push_back(t_segment(align_beg, align_end));

		if (!(mask_region & eWHOLE))
			vec_loci.erase(vec_loci.begin() + 2);
		if (!(mask_region & eINVERSE))
			vec_loci.erase(vec_loci.begin() + 1);
		if (!(mask_region & eMASKED))
			vec_loci.erase(vec_loci.begin() + 0);
		assert(vec_loci.size());


	//____________________________________________________________________________________
	// get strings comprising concatenated domain/free regions
	//____________________________________________________________________________________
		vector<string> proteins;
		for (unsigned j = 0; j < vec_loci.size(); ++j)
			extract_seq_pair_with_loci(proteins, prot1, prot2, vec_loci[j]);

	//____________________________________________________________________________________
	// print gene_id
	//____________________________________________________________________________________
		cout	<< orthologs[i].gene_id1		<< "\t"
				<< orthologs[i].gene_id2;


	//____________________________________________________________________________________
	// print protein lengths
	//____________________________________________________________________________________
		for (unsigned j = 0; j < proteins.size(); j += 2)
			if (proteins[j].length())
				cout << "\t" << proteins[j].length();
			else
				cout << "\t\\N";

	//____________________________________________________________________________________
	// get percentage ids
	//____________________________________________________________________________________
		get_all_regions_percentage_id(proteins);



	//____________________________________________________________________________________
	// get Ka Ks
	//____________________________________________________________________________________
		vector<string> cDNAs;
		for (unsigned j = 0; j < vec_loci.size(); ++j)
			extract_cdna_pair_with_loci(cDNAs, cdna1, cdna2, vec_loci[j]);
		if (!get_all_regions_ka_ks(	cDNAs))
			return 1;


	//____________________________________________________________________________________
	// get bitscore
	//____________________________________________________________________________________
		get_all_regions_bitscores(	prot1, prot2, vec_loci);

	//____________________________________________________________________________________
	// get mouse_to_human_matrix
	//____________________________________________________________________________________
		get_all_regions_matrix(	proteins, freq_lookup_table);

		cout << "\t" << orthologs[i].description;

		cout << "\n";
	}
}











//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	// don't use C streams
	std::ios::sync_with_stdio(false);


	// get command line arguments

	// if no arguments print usage and exit
	if (argc == 1)
	{
		char * my_argv[2] ={ argv[0], "-h"};
		int argc = 2;
		t_program_args args;
		if (!process_command_line_options(argc, my_argv, args))
			return 1;
	}
	t_program_args args;
	if (!process_command_line_options(argc, argv, args))
		return 1;


	my_path = GetPathPart(argv[0]);
	paml_path = args.paml;
	tmp_file_path = my_path + "/seq.tmp";
	yn00_exe_path = paml_path + " " + tmp_file_path;



	short freq_lookup_table[256][256];
	initialize_freq_lookup_table(freq_lookup_table);
	initialize_bitscore_lookup_table(lookup_table, gap_penalty);

	//____________________________________________________________________________________
	// read ortholog data from STDIN
	//____________________________________________________________________________________
	std::vector<ortholog_region> orthologs;
	get_gene_ids(orthologs);

	//____________________________________________________________________________________
	// analyse each ortholog, retrieving the sequence from the specified database
	//____________________________________________________________________________________
	string sql_str = "SELECT sequence1, sequence2, cdna1, cdna2, "
						"gaps1, gaps2, align_start, align_len FROM " +
						args.table;
	if (!analyze_each_ortholog(orthologs, sql_str, args.region, freq_lookup_table,
								args.pg_db.get()))
		return 1;



	return 0;
}

